Two-dimensional Discrete Element Method Research Articles

Assessment of shear band evolution using discrete element modelling

PurposeThe objective of this paper is to quantitatively assess shear band evolution by using two-dimensional discrete element method (DEM).Design/methodology/approachThe DEM model was first calibrated by retrospectively modelling existing triaxial tests. A series of DEM analyses was then conducted with the focus on the particle rotation during loading. An approach based on particle rotation was developed to precisely identify the shear band region from the surrounding. In this approach, a threshold rotation angle ω0 was defined to distinguish the potential particles inside and outside the shear band and an index g(ω0) was introduced to assess the discrepancy between the rotation response inside and outside shear band. The most distinct shear band region can be determined by the ω0 corresponding to the peak g(ω0). By using the proposed approach, the shear band development of two computational cases with different typical localised failure patterns were successfully examined by quantitatively measuring the inclination angle and thickness of shear band, as well as the microscopic quantities.FindingsThe results show that the shear band formation is stress-dependent, transiting from conjugated double shear bands to single shear band with confining stress increasing. The shear band evolution of two typical localised failure modes exhibits opposite trends with increasing strain level, both in inclination angle and thickness. Shear band featured a larger volumetric dilatancy and a lower coordination number than the surrounding. The shear band also significantly disturbs the induced anisotropy of soil.Originality/valueThis paper proposed an approach to quantitatively assess shear band evolution based on the result of two-dimensional DEM modelling.

Study on generation and evolution of base voids in heavy-haul railway tunnels using particle flow method

Most existing freight railway tunnels have experienced basement cracking, damage, and even mud-pumping after a certain period of operation. For heavy-haul railway tunnels with high traffic densities, large capacities, and large axle loads, this situation will be more serious. This paper establishes a simplified model of a single-track heavy-haul railway tunnel using the two-dimensional particle discrete element method. The meso-evolution process of the concrete structure and surrounding rock void at the bottom of the tunnel is studied in the later stages of repeated loading by heavy-haul trains. The difference in the base void between waterless and water-saturated bases is compared. The cracks in the surrounding rock are first generated at the lowest position of the invert center, and run from the middle to the sides and from the top to the bottom. The terminal cracks return to the free surface, exhibiting a parabolic shape that opens upward. There are more cracks in the middle and fewer at the sides. The asynchronous deformation of the surrounding rock and the inverted arch is positively correlated with the number of cracks. The voiding process under water-saturated conditions is similar to that without water, but the depth and width over which the rock mass is crushed and the range of voiding increase exponentially in the water-saturated case, and the void appears earlier. It is anticipated that these research results will contribute to the structural design and operational safety of heavy-haul railway tunnels.

Discrete element modelling of the size-dependent mechanical behavior of Al2O3 ceramic under quasi-static uniaxial compression

The size-dependent mechanical behavior of Al2O3 ceramic under quasit-static uniaxial compression was simulated via two-dimensional discrete element method. The numerical sample of Al2O3 ceramic was established by adopting the linear parallel bond contact model, and the inherent defects was further generated by randomly removing a subset of elementary discs. The uniaxial compressive strength of Al2O3 ceramic samples with different sizes follows the Weibull distribution. With an increase in sample size, the magnitude of the uniaxial compressive strength decreases, while its variability increases. Additionally, with sample size increases, both the Young's modulus and Poisson's ratio decrease. Meanwhile, the variability in Young's modulus and Poisson's ratio is more pronounced in larger samples than in smaller ones. The progressive fracture process of Al2O3 ceramic sample can be categorized into four typical stages, namely linear-elastic stage, crack initiation stage, crack propagation stage and complete failure stage. The final fragment size distributions of Al2O3 ceramic samples of different sizes can be accurately described by the fractal model. Notably, smaller samples tend to exhibit a greater breakage extent compared to larger samples. Finally, a size-dependent statistical equation of uniaxial compressive strength of Al2O3 ceramic is proposed based on the numerical results and the inverse function of the Weibull distribution, which can provide a theoretical basis for the design and production of ceramic materials.

Micro-damage characteristics of cold recycled mixture under freeze–thaw cycles based on discrete-element modeling

Cold recycled mixture with bitumen emulsion (CRME) is an energy-saving and environmentally friendly material. However, the micro-scale damage characteristics under the coupling of water temperature are still unknown and need to be investigated. In this paper, based on the X-ray CT section, the two-dimensional discrete element method (DEM) was used to precise reconstruct the CRME under different freeze–thaw (F-T) cycles. Simultaneously, the fracture behavior, strength, and structural damage characteristics of CRME under F-T cycles were studied by DEM and digital image correlation technology. The results show that F-T cycle significantly damage the formation of strength force chain in CRME. After undergoing 20F-T cycles, the specimen exhibited tensile failure without any obvious crack penetration. The cracks required of the specimen after the F-T cycles to reach the peak load were less than the non-F-T. The stress–strain concentrations were reinforced, making the specimen more susceptible to unloading during the later stages of testing after F-T cycles. The repeated F-T cycles damaged the overall microstructure of the mixture, resulting in less noticeable tensile effects during the splitting test. Utilizing a precise DEM modeling method proved to be practical and effective in investigating the micro-damage characteristics of cold recycled mixture.

DEM study on the effect of particle shape on the shear behaviour of granular materials

This study investigates the effects of particle convexity, sphericity and aspect ratio (AR) on the behaviour of sheared granular materials using two-dimensional discrete element method simulations. Isotropic, dense and loose assemblies with different particle shapes were prepared and subjected to drained shearing via biaxial compression until the critical state was reached. Macroscopic characteristics such as strength and dilatancy are presented. The factors underlying the macroscopic behaviour are then investigated by considering the coordination number, fabric anisotropy, particle moment, friction mobilisation at contacts and particle rotation. For the range of shapes considered here, the data indicate that the shear strength decreases as particle convexity and sphericity increases while the shear strength increases with increasing AR. The shear strength and convexity are weakly correlated, however a stronger correlation is observed between AR and strength. The volumetric strain at large strains tends to increase with increasing AR. There is a stronger correlation between the critical state strength and both the critical state coordination number and the critical state mechanical void ratio than there is between the critical state void ratio and the critical state strength. The contact fabric anisotropy, the magnitude of the moment transmitted by particles and the friction mobilised at the contacts are important factors underlying strength. The critical state strength increases as both the mean particle moment and the mean mobilised friction increased. Analysis of particle rotation provides insights into the response of the granular materials to shearing.

DEM Analysis of Dynamic Evolutions of Lateral Soil Arching in Sandy Soil

For landslide mitigation, stabilizing piles are widely adopted in geotechnical engineering practice. The formation of soil arching between stabilizing piles can greatly mitigate the lateral spreading of sandy soil, and complementarily, evidence demonstrates that the sandy soil properties have a significant influence on the soil arching process. Thus, this paper focuses on investigating the lateral arching evolution behaviors in both dense and loose sands. A discrete element method (DEM) involving the servo-mechanism that plays a role in quantitatively controlling the translational velocity of selected walls when a desired force is to be applied or maintained was carried out to simulate the dynamic evolutions of the lateral soil arching process. The double-arch model is proposed to illustrate the end-bearing arch and frictional arch to analyze soil arching effects separately. Results indicate that the macrobehaviors and microbehaviors of lateral soil arching in dense and loose sands are initially distinct but ultimately similar. The arching mechanism is also illustrated by fabric analysis and coordination numbers. In general, the lateral arching evolutions in both dense and loose sands can be divided into three evolutionary stages as the loading displacement increases, which are similar to the strain-hardening and strain-softening phenomenon, respectively. Similarly, a unique arching process of forming-breaking-forming is identified by the use of particle trace. Moreover, multiple two-dimensional DEM models at different depths are superimposed vertically to study the three-dimensional soil arching. Based on the parametric studies, the loading velocity and pile spacing ratio have a preponderant influence on the development of soil arching, whereas the influence of the depth and confining stress is almost negligible. Additionally, the distribution of horizontal thrust has been determined by monitoring the stabilizing pile-based arching foot at different locations and the distribution of shear resistance is found to be approximately triangular.

The Sensitivity of Micro—Macro Mechanical Behaviour of Sand to the Inter-Particle Properties

Sand is a particulate material but is treated as a continuum solid in some engineering analyses. This approach is proven to be acceptable when dealing with geotechnical structures, provided an adequate factor of safety is applied so that there is no risk of failure. However, the continuum approach does not account for the effect of interparticle forces on the micro–macro behaviour of sand. Sand could be modelled as a particulate material using the discrete element method (DEM), taking into account its discrete nature. This paper shows how the microscopic contact properties between the idealised sand particles influence the macro-mechanical behaviour, highlighting the development of the fabric as the soil approaches failure. Thirty DEM biaxial tests were performed to study the sensitivity of the macro–micro mechanical properties of sand to the inter-particle properties of an idealised sand particle. The conditions of these simulations were the same (e.g., particle size distribution, number of particles, porosity after radius enlargement, boundary conditions, and rate of loading). The sensitivity of the pre-peak, peak, and post-peak behaviour of these simulations to the inter-particle properties of an idealised sand particle was studied. Two extra DEM biaxial tests under different confining pressures were performed to verify the cohesionless nature of the synthetic material used for this study. Since a two-dimensional DEM is used for this study, a detailed approach to interpret the results assuming either a plane strain or a plane stress situation was discussed. This study highlighted the critical inter-particle properties and the range over which these influence macro-mechanical behaviour. The results show that Young’s modulus is mainly dependent on the normal contact stiffness, and peak stress and the angle of internal friction are greatly dependent on the inter-particle coefficient of friction, while Poisson’s ratio and volumetric behaviour of particulate sand are dictated mainly by shear contact stiffness. A set of relationships were established between inter-particle properties and macro-machinal parameters such as Young’s modulus, Poisson’s ratio, and angle of internal friction. The elastoplastic parameters obtained from these tests are qualitatively in agreement with the typical medium and dense sand behaviour.

Study of Strain-rate Effect in Two-dimensional Biaxial Test on Granular Material using Discrete Element Method

Numerical study to investigate the effect of various strain rates on global friction angle in the sand has been performed. Granular material behavior is influenced by several factors, among others: pack configuration, grain macro and micro roughness, confinement pressure, loading rate, etc. Sand is a granular material composed of discrete particles that the most refined microscopic techniques are needed to study its mechanical properties. In Indonesia, research related to the Discrete Element Method is still very limited. The two-dimensional discrete element method is capable to calculate the motion and interparticle contacts of large number of small particles, and each particle is modeled as a rigid circular element. The study started with the validation of the DEM model using YADE. Particles with a local friction angle of 35° are arranged in a closed rectangular box (frictionless wall). The number of particles used in this model validation simulation is 1000 sphere-type particles with monodisperse particle gradations. The simulation was done by a drained biaxial test with confinement stress of 100 kPa, thereafter varying strain rate are applied. Based on the deviatoric stress–axial strain curve from YADE, the result can be depicted on the Mohr circle to obtain the value of the global friction angle. It is found that the different value of strain rate affects its global friction angle. Increasing the value of the strain rate can increase the material global friction angle, which increases the strain rate from 1% to 5%, 10%, 25%, and 50% will increase its global friction angle by 5%, 5%, 14%, and 18%, respectively.

Understanding of penetrating resistance on an inclined plate to cohesive soil using two-dimensional discrete element method

In order to know the effect of varied grousers on the reaction of soil, the reaction of soil by a plate at different angles of inclination using the two-dimensional discrete element method (DEM) has been investigated. Either the simulation or the experiment, a flat iron plate with 5 × 77 × 100 (mm) was used, and the inclination angle between the plate and the soil surface ranged from 0 to 60 degrees with 15-degree intervals. According to the result, the cohesive soil model introduced in the 2D DEM could successfully replicate the experimental results considering the effects of different inclination angles of the plate and the values of the root mean square error between the simulation and the experiment were found as 6.62, 36.60, 58.63, 60.84, and 66.25 N for inclination angles from 0 to 60 degrees in the horizontal direction. Moreover, the RMSE values of 127.91, 216.06, 75.88, 110.40, and 25.49 N were obtained for the vertical direction. Overall, a minimal discrepancy was observed between the results of DEM simulation and those of experiments in the vertical component of reaction force when the inclination angle of the plate was increased.

Effects of confining pressure and cutting sequence on the cobalt-rich crust cutting mechanism by using conical picks in ocean environments

Confining pressure and cutting sequence of conical picks directly affect the ocean mining efficiency. The simulation model of cobalt-rich crust cutting with two-dimensional discrete element method is established to simulate the whole process of crack growth. According to the obtained crack growth patterns and the rules of fracture-block formation, four typical fracture modes are summarized, and the fracture mechanism is further studied from three aspects: the efficiency of cobalt-rich crust cutting, the ability of crack growth and the size of fracture-block. The results show that with the increase of confining pressure, the cutting efficiency and crack growth ability decrease, and the cutting efficiency of simultaneous loading is larger than that of the sequential loading, while the crack growth ability is just the opposite. In addition, the size of fracture-block is related to fracture modes.

On the External Failure Surface in PSEs: Numerical and Theoretical Methods

The evolution of soil arching, involving internal and external failure surfaces, is of significance to the load transfer mechanism of pile-supported embankments (PSEs), and external failure surfaces are generally observed in cases of greater embankment height. In this study, the evolution of external failure surfaces was investigated by both numerical and theoretical methods. To begin with, numerical simulations of trapdoor tests were carried out by a two-dimensional discrete element method. The influences of two key parameters (i.e., embankment height and net spacing between piles) on the development of external failure surfaces were emphasized. The measured coefficient of lateral earth pressure around the external failure surface was close to the predicted coefficient of active earth pressure. Then, a theoretical solution considering inclined external failure surfaces, more realistic compared to those adopted in Terzaghi’s method, was proposed. Compared with Terzaghi’s method, the proposed solution exhibited better consistency with laboratory observations, especially when external failure surfaces were significantly observed.

Loading rate dependence of granular cushion under rockfall impact and proposal of an analytical model for impact force estimation

The cushioning behaviour of granular materials, such as sand and gravel, has been studied for the rational design and performance evaluation of a rock shed intended for rockfall protection. The shape of the collision surface of a falling mass strongly affects the impact force waveform of the falling mass. However, the control factors for generating impact force waveform have not been fully clarified. Previous studies have shown that the maximum transmitted impact force acting on the top plate of a rock shed through a cushioning material depends not only on the maximum impact force acting on the falling mass but also on its duration. To accurately predict the maximum transmitted impact force acting on a structure, understanding the generation mechanism of the impact force waveform acting on the falling mass is necessary. These control factors were investigated through the two-dimensional (2D) discrete element method (DEM) analysis. It was shown that the mean principal stress of particles around a falling mass increased with the loading velocity of the falling mass. Based on this result, an analytical model for estimation the impact force waveform of the falling mass was proposed. The validity of the proposed method and the control factors estimated via simulations was confirmed by comparing the results obtained by the proposed method with 2D DEM simulation results and experimental results. This study contributes to the elucidation of the deformation mechanism of soil under high-speed loading and development of rational countermeasures considering the deformation behaviour of a sand cushion.

Effects of fibroblast on electromechanical dynamics of human atrial tissue-insights from a 2D discrete element model.

Roughly 75% of normal myocardial tissue volume is comprised of myocytes, however, fibroblasts by number are the most predominant cells in cardiac tissue. Previous studies have shown distinctive differences in cellular electrophysiology and excitability between myocytes and fibroblasts. However, it is still unclear how the electrical coupling between the two and the increased population of fibroblasts affects the electromechanical dynamics of cardiac tissue. This paper focuses on investigating effects of fibroblast-myocyte electrical coupling (FMEC) and fibroblast population on atrial electrical conduction and mechanical contractility by using a two-dimensional Discrete Element Method (DEM) model of cardiac tissue that is different to finite element method (FEM). In the model, the electro-mechanics of atrial cells are modelled by a biophysically detailed model for atrial electrical action potentials and myofilament kinetics, and the atrial fibroblasts are modelled by an active model that considers four active membrane ionic channel currents. Our simulation results show that the FMEC impairs myocytes’ electrical action potential and mechanical contractibility, manifested by reduced upstroke velocity, amplitude and duration of action potentials, as well as cell length shortening. At the tissue level, the FMEC slows down the conduction of excitation waves, and reduces strain of the tissue produced during a contraction course. These findings provide new insights into understandings of how FMEC impairs cardiac electrical and mechanical dynamics of the heart.

Lethe-DEM: an open-source parallel discrete element solver with load balancing

Approximately $${75}\%$$ of the raw material and $${50}\%$$ of the products in the chemical industry are granular materials. The discrete element method (DEM) provides detailed insights of phenomena at particle scale, and it is therefore often used for modeling granular materials. However, because DEM tracks the motion and contact of individual particles separately, its computational cost increases nonlinearly $$O(n_\mathrm{p}\log (n_\mathrm{p}))$$ – $$O(n_\mathrm{p}^2)$$ (depending on the algorithm) with the number of particles ( $$n_\mathrm{p}$$ ). In this article, we introduce a new open-source parallel DEM software with load balancing: Lethe-DEM. Lethe-DEM, a module of Lethe, consists of solvers for two-dimensional and three-dimensional DEM simulations. Load balancing allows Lethe-DEM to significantly increase the parallel efficiency by $$\approx {25}$$ – $${70}\%$$ depending on the granular simulation. We explain the fundamental modules of Lethe-DEM, its software architecture, and the governing equations. Furthermore, we verify Lethe-DEM with several tests including analytical solutions and comparison with other software. Comparisons with experiments in a flat-bottomed silo, wedge-shaped silo, and rotating drum validate Lethe-DEM. We investigate the strong and weak scaling of Lethe-DEM with $${1}\le n_\mathrm{c} \le {192}$$ and $${32}\le n_\mathrm{c} \le {320}$$ processes, respectively, with and without load balancing. The strong-scaling analysis is performed on the wedge-shaped silo and rotating drum simulations, while for the weak-scaling analysis, we use a dam-break simulation. The best scalability of Lethe-DEM is obtained in the range of $${5000}\le n_\mathrm{p}/n_\mathrm{c} \le {15{,}000}$$ . Finally, we demonstrate that large-scale simulations can be carried out with Lethe-DEM using the simulation of a three-dimensional cylindrical silo with $$n_\mathrm{p}={4.3}\times 10^6$$ on 320 cores.

Particle impact on a cohesive granular media.

We investigate numerically the impact process of a particle of diameter d and velocity V_{i} onto a cohesive granular packing made of similar particles via two-dimensional discrete element method simulations. The cohesion is ensured by liquid bridges between neighboring particles and described by short range attraction force based on capillary modeling. The outcome of the impact is analyzed through the production of ejected particles from the packing, referred to as the splash process. We quantify this production as a function of the impact velocity for various capillary strength Γ and liquid content Ω. The numerical data indicate that the splash process is modified when the dimensionless cohesion number Co=6Γ/ρ_{p}gd^{2} (where ρ_{p} is the particle density, d its diameter, and g the gravitational acceleration) exceeds a critical value of the order of the unity. Above this value, we highlight that the ejection process is triggered above a threshold impact Froude number, Fr=V_{i}/sqrt[gd], which depends both on Γ and Ω and scales as Γ^{β}Ω^{δ}, where the values of the exponents are found close to 1/2 and 1/6, respectively, and can be derived from rational physical arguments. Importantly, we show that, above the threshold, the number of splashed particles follows a linear law with the impact Froude number as in the cohesionless case.

Sheared granular matter and the empirical relations of seismicity.

The frictional instability associated with earthquake initiation and earthquake dynamics is believed to be mainly controlled by the dynamics of fragmented rocks within the fault gauge. Principal features of the emerging seismicity (e.g., intermittent dynamics and broad time and/or energy scales) have been replicated by simple experimental setups, which involve a slowly driven slider on top of granular matter, for example. Yet these setups are often physically limited and might not allow one to determine the underlying nature of specific features and, hence, the universality and generality of the experimental observations. Here, we address this challenge by a numerical study of a spring-slider experiment based on two-dimensional discrete element method simulations, which allows us to control the properties of the granular matter and of the surface of the slider, for example. Upon quasistatic loading, stick-slip-type behavior emerges which is contrasted by a stable sliding regime at finite driving rates, in agreement with experimental observations. Across large parameter ranges for damping, interparticle friction, particle polydispersity, etc., the earthquake-like dynamics associated with the former regime results in several robust scale-free statistical features also observed in experiments. At first sight, these closely resemble the main empirical relations of tectonic seismicity at geological scales. This includes the Gutenberg-Richter distribution of event sizes, the Omori-Utsu-type decay of aftershock rates, as well as the aftershock productivity relation and broad recurrence time distributions. Yet, we show that the correlations associated with tectonic aftershocks are absent such that the origin of the Omori-Utsu relation, the aftershock productivity relation, and Båth's relation in the simulations is fundamentally different from the case of tectonic seismicity. This, we believe, is mainly due to a lack of macroscale relaxation processes that are closely tied to the generation of real aftershocks. We argue that the same is true for previous laboratory experiments.

DEM simulation of rock creep in tunnels using Rate Process Theory

The time-dependent (creep) behaviour of rocks affects the safety and stability of tunnels excavated in weak rocks and at great depths. Several theories have been proposed to simulate the creep deformation in rock; i.e., the progressive time-dependent damage that rocks (or other materials) exhibit under constant stress. However, most of these theories do not capture the accelerating strains associated to tertiary creep and leading to rock failure. In this research, the Rate Process Theory (RPT), combined with the Discrete Element Method (DEM), are used to simulate rock creep deformation in deep tunnels. To do that, two-dimensional (2D) DEM tunnel models are built using particles, with their interactions being simulated by a hybrid mixture of the linear and flat joint contact models. The RPT is incorporated into such models by a user-defined Visual C++ script that modifies their friction coefficients during the DEM simulation, depending on the relative velocity between particles. Numerical results show, for the first time, that the joint RPT-DEM approach is able to reproduce all stages of tunnel convergences associated to rock creep, including tertiary creep; and that it can reproduce the rock damage associated to such creep strains.

Effect of relative density of 2D granular materials on the arching effect through numerical trapdoor tests

The relative density and the fill height of granular materials influence the arching effect significantly, from either the overall behavior of granular materials or the formation of the arching effect. In this study, a comprehensive comparison of the arching effect at four fill heights with three relative densities is conducted by numerical trapdoor tests using the two-dimensional (2D) discrete element method. The results indicate that the stick–slip behavior of the arching ratio (the ratio of the average vertical stress on the trapdoor to the overburden stress; the ratio decreases to a minimum value, followed by recovering to the ultimate value) is more obvious for the case with lower fill height at a given relative density or that with higher density at a given fill height. The partial and the full arching can be distinguished by whether the loosen zone propagates to the surface. Two different critical states are identified for the partial and the full arching, respectively. These critical states are insensitive to the relative density. A characteristic height is identified around twice the width of the trapdoor, at which the critical state of the partial arching transfers to that of the full arching with the increase of the relative density.

DEM numerical investigation on the effect of back rake angle of PDC bit in drilling hard rock strata

Polycrystalline diamond compact (PDC) drill bits used in the Deep Borehole Placement (DBP) approach for nuclear waste disposal are subjected to highly competent, low-porosity, and low-permeability rock mass that cause severe wear and damage. To augment the penetration rate and minimize the wear-induced cost, a novel interaction rock cutting model with breakable granite and unbreakable cutter was established using two-dimensional discrete element method (2D-DEM). The particle element was used to simulate the granite sample, and the micro parameters of the contact elements between them were calibrated based on the uniaxial compression and Brazilian tests. While the rigid block element was used to reproduce the single PDC cutter including its elastic properties and geometric structure (i.e., back rake angle). The cutting process was simulated by applying normal force and tangential velocity on the single cutter to explore the effect of the back rake angle (5, 10, 15, 20 and 25°) on the stress condition of the cutter and penetration rate. The simulation results show that with the increase of the back rake angle, the maximum shear stress magnitude and stress concentration area decrease simultaneously, in according with the decrease of material removal rate, which do good help to decrease the wear and unexpected damage of bit cutter. In addition, the cutting process can be divided into the effective cutting stage when the elastic energy continuously stores and the fragmentation initiation stage when the stored energy suddenly releases and numerous cracks generate. In long term, the established DEM model can be used to optimize the structure of the PDC bit.

Discrete Element Modeling of Vibration Compaction Effect of the Vibratory Roller in Roundtrips on Gravels

Abstract This paper aims to study the vibration compaction mechanism of the vibratory roller on gravels using a two-dimensional discrete element method. The roadbed model was established by gravel particles with irregular shapes, which was closer to reality. The performance parameters of the vibratory roller, such as operating frequency and rolling velocity, were investigated to explore their influences on the operating efficiency of the vibratory roller in roundtrips. The frequencies of 15 Hz and 17 Hz were proved to be the optimal frequency and resonance frequency in the current simulations, respectively. The vibratory roller could achieve a better vibration compaction effect with less power consumption at the optimal frequency. In addition, the number of roundtrips and power consumption should be considered in the selection of the optimal rolling velocity. The movement direction and the contact force distribution of gravels were illustrated by the displacement field, velocity field, as well as the contact force chains. Our results provide a better understanding of the mechanical behavior of gravel particles and their interactions with the vibratory roller.